home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / Lyapunov / lyap / lyap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-07-31  |  11.6 KB  |  571 lines
  1. /*
  2.  * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
  3.  */
  4.  
  5. /* This program is Copyright (C) 1991, Garrett A. Wollman.  This
  6.  * program may be modified and distributed for any purpose and without
  7.  * fee, provided that this notice remains on all such copies
  8.  * unaltered.  Binary distributions not including this source file are
  9.  * prohibited.  Modified versions must be marked with the name of the
  10.  * modifier and the date of modification.
  11.  *
  12.  * Under no circumstances shall Garrett A. Wollman or the University
  13.  * of Vermont and State Agricultural College be held liable for any
  14.  * damages resulting from the use or misuse of this program, whether
  15.  * the author is aware of such possibility or not.  This program is
  16.  * warranted solely to occupy disk space.
  17.  *
  18.  * I'm sorry for all this legal garbage, but it is sadly necessary
  19.  * these days.
  20.  */
  21.  
  22. #include <math.h>
  23. #include <time.h>
  24. #ifndef SUN_BROKEN_STDLIB
  25. #include <stdlib.h>
  26. #endif
  27. #include <stdio.h>
  28. #ifdef sgi
  29. #ifndef NOMULTI
  30. #include <sys/types.h>
  31. #include <sys/prctl.h>
  32. #include <sys/wait.h>
  33.  
  34. #define MULTIPROC
  35. #define DO_SPROC do_sproc
  36. #define MAXPROCS 8
  37. #endif
  38. #endif
  39.  
  40. /* other multiprocessing hosts please add your info here */
  41.  
  42.  
  43. /*
  44.  * Global variables
  45.  */
  46. #ifdef MULTIPROC
  47. int nprocs = 1;
  48. #endif
  49.  
  50. char *whoami;
  51.  
  52. extern int getopt(int,char **,const char *);
  53. extern char *optarg;
  54. extern int optind;
  55. extern int opterr;
  56. extern void perror(const char *);
  57.  
  58. #ifdef SUN_BROKEN_STDLIB
  59. extern volatile void exit(int);
  60. extern void *malloc(int);
  61. extern int atoi(const char *);
  62. #endif
  63.  
  64. extern int strlen(const char *);
  65.  
  66. #ifdef NO_STDIO_PROTOS
  67. extern int fwrite(const char *,int,int,FILE *);
  68. extern int fputc(int,FILE *);
  69. extern int fputs(const char *,FILE *);
  70. extern int pclose(FILE *);
  71. extern void fclose(FILE *);
  72. extern int fprintf(FILE *,const char *,...);
  73. extern int sscanf(char *,const char *,...);
  74. #endif
  75.  
  76. #ifdef MULTIPROC
  77. #define OPTIONS "r:c:v:m:d:s:g:5pt#:o:l:"
  78. #else
  79. #define OPTIONS "r:c:v:m:d:s:g:5pto:l:"
  80. #endif
  81.  
  82. int rows = 512,cols = 512;
  83. double amin,amax;
  84. double bmin,bmax;
  85. double aincr,curA;
  86. double bincr;
  87.  
  88. #define nColors 256        /* don't even try to change this */
  89. struct {
  90.   unsigned char red;
  91.   unsigned char green;
  92.   unsigned char blue;
  93. } colormap[nColors];
  94.  
  95. int Settle = 600;
  96. int Dwell = 2000;
  97. int *Rvec;
  98. double minLya = -5;
  99. int showpos = 0;
  100. double initGuess = 0.5;
  101.  
  102. /*
  103.  * Colormap functions
  104.  */
  105.  
  106. /*
  107.  * set up our initial shades of grey
  108.  */
  109. void init_colormap(void) {
  110.   int i;
  111.   for(i=0; i < nColors; i++)
  112.     colormap[i].red = colormap[i].green = colormap[i].blue = i;
  113. }
  114.  
  115. void read_colormap(const char *name) {
  116.   FILE *cmap;
  117.   char buf[255];        /* magic number */
  118.   int a,b,c,i;
  119.   
  120.   cmap = fopen(name,"r");
  121.   if(!cmap) {
  122.     perror(name);
  123.     exit(-1);
  124.   }
  125.  
  126.   i = 0;
  127.   while(!feof(cmap) && i < nColors) {
  128.     fgets(buf,sizeof buf,cmap);
  129.     /*
  130.      * unparseable lines are comments
  131.      * as is anything after the <r g b> value on a single line
  132.      *
  133.      * per Fractint 15.1 standard
  134.      */
  135.     if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
  136.       continue;
  137.  
  138.     colormap[i].red = a;
  139.     colormap[i].green = b;
  140.     colormap[i].blue = c;
  141.     i++;
  142.   }
  143.   fclose(cmap);
  144. }
  145.  
  146.  
  147. /*
  148.  * Lyapunov exponent and coloring calculation
  149.  *
  150.  * original function by Ed Kubaitis, used by permission
  151.  */
  152.  
  153. /*
  154.  * Speedup...
  155.  *
  156.  * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
  157.  * is, by a principle from high-school algebra, the same as
  158.  * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
  159.  * by unrolling the dwell loop by four (hence the rounding in main())
  160.  * and accumulating by multiplication before taking the log.  Since
  161.  * log() is extremely expensive, this should lead to a considerable
  162.  * speedup, while still allowing complete flexibility in selecting
  163.  * dwell values.
  164.  *
  165.  * Presumably, one might unroll this loop even further, with a 
  166.  * concomitant increase in textual complexity.
  167.  */
  168.  
  169. int lya(double a,double b) {
  170.   double t = 0;
  171.   double r, lya;
  172.   int d;
  173.   double acc;
  174.   double x;
  175.  
  176.   x = initGuess;        /* initialize x */
  177.  
  178.   for(d=0; d < Settle; d++) {
  179.     r = (Rvec[d]) ? b : a;
  180.     x = r*x*(1-x);
  181.   }
  182.  
  183. #if 1
  184.   for(d=0; d <= Dwell; d+= 4) {
  185.     r = (Rvec[d]) ? b : a;
  186.     x = r * x * (1-x);
  187.     acc = fabs(r - 2*r*x);
  188.  
  189.     r = (Rvec[d+1]) ? b : a;
  190.     x = r * x * (1-x);
  191.     acc *= fabs(r - 2*r*x);
  192.  
  193.     r = (Rvec[d+2]) ? b : a;
  194.     x = r * x * (1-x);
  195.     acc *= fabs(r - 2*r*x);
  196.  
  197.     r = (Rvec[d+3]) ? b : a;
  198.     x = r * x * (1-x);
  199.     t += log(acc * fabs(r - 2*r*x));
  200.  
  201.     if(fabs(x-0.5) < 1e-20) {
  202.       d += 3;
  203.       break;
  204.     }
  205.   }
  206. #else
  207.   for(d = 0; d <= Dwell; d += 4) {
  208.     r = (Rvec[d]) ? b : a;
  209.     x = r * x * (1-x);
  210.     t += log(fabs(r-2*r*x));
  211.  
  212.     r = (Rvec[d+1]) ? b : a;
  213.     x = r * x * (1-x);
  214.     t += log(fabs(r-2*r*x));
  215.  
  216.     r = (Rvec[d+2]) ? b : a;
  217.     x = r * x * (1-x);
  218.     t += log(fabs(r-2*r*x));
  219.  
  220.     r = (Rvec[d+3]) ? b : a;
  221.     x = r * x * (1-x);
  222.     t += log(fabs(r-2*r*x));
  223.  
  224.     if(fabs(x-0.5) < 1e-20) {
  225.       d += 3;
  226.       break;
  227.     }
  228.   }
  229. #endif
  230.  
  231.   lya = (t * 1.442695041)/d;
  232.  
  233.   if(showpos)            
  234.     lya *= -1;
  235.   if(lya < minLya)        /* sanity check! */
  236.     lya = minLya;
  237.  
  238.   if(lya < 0) {
  239.     return (int)(lya / minLya * (double)nColors);
  240.   } else {
  241.     return 0;
  242.   }
  243. }
  244.  
  245.  
  246. /*
  247.  * Usage message
  248.  */
  249. #ifdef __GNUC__
  250. volatile
  251. #endif
  252. void usage(void) {
  253.   fprintf(stderr,
  254.       "%s: usage:\n\n"
  255.       "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
  256. #ifdef MULTIPROC
  257.       "[-# processors] "
  258. #endif
  259.       "[-o file] [-l colormap]\n\t"
  260.       "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
  261.       "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  262.   exit(-1);
  263. }
  264.  
  265.  
  266. /*
  267.  * SGI-specific multiprocessing code
  268.  */
  269. #ifdef sgi
  270. #ifndef NOMULTI
  271. struct info {
  272.   int *row;
  273.   double curA;
  274. };
  275.  
  276. void real_doit(double curA,int *row) {
  277.   int onCol;
  278.   double curB;
  279.  
  280.   for(onCol = cols - 1,curB = bmax; onCol; onCol--,curB -= bincr) {
  281.     row[onCol] = lya(curA,curB);
  282.   }
  283. }
  284.  
  285. void doit(void *info) {
  286.   struct info *myinfo = info;
  287.  
  288.   real_doit(myinfo->curA,myinfo->row);
  289.   exit(0);
  290. }
  291.  
  292. void do_sproc(int *manyRows[]) {
  293.   int pids[MAXPROCS];
  294.   struct info infos[MAXPROCS];
  295.   int i;
  296.   double aval = curA;
  297.  
  298.   for(i=0; i < nprocs - 1; i++,aval -= aincr) {
  299.     infos[i].row = manyRows[i];
  300.     infos[i].curA = aval;
  301.     pids[i] = sproc(doit,PR_SADDR,(void *)&infos[i]);
  302.   }
  303.   real_doit(aval,manyRows[i]);
  304.  
  305.   /*
  306.    * We know we started (nprocs - 1) kids, so we wait for (nprocs - 1) kids
  307.    * to die.
  308.    */
  309.   for(i=0; i < nprocs - 1; i++) {
  310.     (void)wait((int *)0);
  311.   }
  312. }
  313. #endif
  314. #endif
  315.  
  316.  
  317. int main(int argc,char **argv) {
  318.   int c,i,j;
  319.   char *viewprogram = NULL;
  320.   char *outname = NULL;
  321.   FILE *outfile;
  322.  
  323.   char *bits;
  324.  
  325.   double curB;
  326.   int onRow,onCol;
  327. #ifdef MULTIPROC
  328.   int *manyRows[MAXPROCS];
  329. #endif
  330.   int *oneRow;
  331.  
  332.   int usepgm = 0;
  333.   int outtext = 0;
  334.  
  335.   init_colormap();
  336.   opterr = 0;
  337.  
  338.   whoami = *argv;
  339.   while((c = getopt(argc,argv,OPTIONS)) != EOF) {
  340.     switch(c) {
  341.     case 'r':
  342.       rows = atoi(optarg);
  343.       if(!rows)
  344.     rows = 512;
  345. /*
  346.     usage();
  347. */
  348.       break;
  349.  
  350.     case 'c':
  351.       cols = atoi(optarg);
  352.       if(!cols)
  353.     cols = 512;
  354. /*
  355.     usage();
  356. */
  357.       break;
  358.  
  359. #ifdef MULTIPROC
  360.     case '#':
  361.       nprocs = atoi(optarg);
  362.       if(!nprocs)
  363.     usage();
  364.       break;
  365. #endif
  366.  
  367.     case 'o':
  368.       outname = optarg;
  369.       viewprogram = NULL;
  370.       break;
  371.  
  372.     case 'v':
  373.       viewprogram = optarg;
  374.       outname = NULL;
  375.       break;
  376.  
  377.     case 'l':
  378.       read_colormap(optarg);
  379.       break;
  380.  
  381.     case 'm':
  382.       if(!sscanf(optarg,"%lf",&minLya))
  383.     usage();
  384.       break;
  385.  
  386.     case 'd':
  387.       Dwell = atoi(optarg);
  388.       if(!Dwell)
  389.     usage();
  390.       /*
  391.        * By forcing the Dwell to be a multiple of four, we can use the
  392.        * speedup noted in lya().  We *could* just do this silently,
  393.        * by if it might make a difference to someone...
  394.        */
  395.       if(Dwell % 4) {
  396.     fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
  397.         Dwell,
  398.         Dwell += 4 - (Dwell % 4));
  399.       }
  400.       break;
  401.  
  402.     case 's':
  403.       Settle = atoi(optarg);
  404.       if(!Settle)        /* for settle close to 0, use 1 */
  405.     usage();
  406.       break;
  407.  
  408.     case 'g':
  409.       if(!sscanf(optarg,"%lf",&initGuess))
  410.     usage();
  411.       break;
  412.  
  413.     case '5':
  414.       usepgm = 1;
  415.       break;
  416.  
  417.     case 'p':
  418.       showpos = 1;
  419.       break;
  420.  
  421.     case 't':
  422.       outtext = 1;
  423.       break;
  424.  
  425.     default:            /* includes '?' (and '#' when !MULTIPROC) */
  426.       usage();
  427.     }
  428.   }
  429.  
  430.   if(!argv[optind])
  431.     usage();
  432.  
  433.   /*
  434.    * Note to language nit-pickers... The code below is legal because
  435.    * the || operator is a sequence point.
  436.    */
  437.   if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
  438.      !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
  439.      !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
  440.      !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
  441.     usage();
  442.  
  443.   bits = argv[optind];
  444.  
  445.   if(Dwell < Settle)
  446.     Dwell = Settle;        /* a bit of sanity */
  447.  
  448.   j = strlen(bits);
  449.  
  450.   Rvec = (int *)malloc((Dwell+1) * sizeof *Rvec);
  451.   if(!Rvec) {
  452.     fputs("Cannot allocate space for Markus vector!\n",stderr);
  453.     exit(-1);
  454.   }
  455.  
  456.   for(i=0; i <= Dwell; i++) {
  457.     if(bits[i % j] == 'a')
  458.       Rvec[i] = 0;
  459.     else if(bits[i % j] == 'b')
  460.       Rvec[i] = 1;
  461.     else {
  462.       fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
  463.           (int)(unsigned char)bits[i % j], bits[i % j]);
  464.       exit(-1);
  465.     }
  466.   }
  467.  
  468. #ifdef MULTIPROC
  469.   for(i=0; i < nprocs; i++) {
  470.     manyRows[i] = malloc(cols * sizeof *manyRows[i]);
  471.     if(!manyRows[i]) {
  472.       fputs("Cannot allocate space for row buffers!\n",stderr);
  473.       exit(-1);
  474.     }
  475.   }
  476. #else 
  477.   oneRow = (int *)malloc(cols * sizeof *oneRow);
  478.   if(!oneRow) {
  479.     fputs("Cannot allocate space for single-row buffer!\n",stderr);
  480.     exit(-1);
  481.   }
  482. #endif
  483.  
  484.   aincr = (amax - amin) / rows;
  485.   bincr = (bmax - bmin) / cols;
  486.  
  487.   if(viewprogram) {
  488.     outfile = popen(viewprogram,"w");
  489.     if(!outfile) {
  490.       fprintf(stderr,"Couldn't popen(\"%s\",\"w\")!\n",viewprogram);
  491.       exit(-1);
  492.     }
  493.   } else if(outname) {
  494.     outfile = fopen(outname,"w");
  495.     if(!outfile) {
  496.       perror(outname);
  497.       exit(-1);
  498.     }
  499.   } else {
  500.     outfile = stdout;
  501.   }
  502.  
  503.   /*
  504.    * Print P[PG]M header...
  505.    */
  506.  /* fprintf(outfile,"P%d %d %d\n",(usepgm?2:3) + (outtext?0:3),cols,rows);
  507.   fprintf(outfile,"%d\n",nColors-1);*/
  508.   putc('\0',outfile);  putc('\0',outfile);
  509.   putc('\0',outfile);  putc('\0',outfile);
  510.   putc('\0',outfile);  putc('\0',outfile);
  511.   putc('\0',outfile);  putc('\0',outfile);
  512.   putc('\0',outfile);  putc('\0',outfile);
  513.   putc('\0',outfile);  putc('\0',outfile);
  514.   fprintf(stderr, "Sizes: %d %d %d %d\n",cols%256,cols/256,rows%256,rows/256);
  515.   putc((unsigned char)(cols%256),outfile);
  516.   putc((unsigned char)(cols/256),outfile);
  517.   putc((unsigned char)(rows%256),outfile);
  518.   putc((unsigned char)(rows/256),outfile);
  519.   putc('\0',outfile);  putc('\0',outfile);
  520.  
  521.  
  522. #ifndef MULTIPROC
  523.   for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
  524.     for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
  525.       oneRow[onCol] = lya(curA,curB);
  526.     }
  527.     
  528. #else
  529.   for(onRow = 0, curA = amax; onRow < rows;
  530.       curA -= aincr*nprocs, onRow += nprocs) {
  531.     DO_SPROC(manyRows);
  532.  
  533.     for(j = 0; j < nprocs; j++) {
  534.       oneRow = manyRows[j];
  535. #endif
  536.     for(i=0; i < cols; i++) {
  537.       unsigned char c;
  538.       c = (unsigned char)oneRow[i];
  539.       if(outtext) {
  540.     if(usepgm)
  541.       fprintf(outfile,"%d ",c);
  542.     else
  543.       fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
  544.           colormap[c].blue);
  545.       } else if(usepgm) {
  546.     fwrite((char *)&c,sizeof c,1,outfile);
  547.       } else {
  548.       putc(colormap[c].red,outfile);
  549.       putc(colormap[c].green,outfile);
  550.       putc(colormap[c].blue,outfile);
  551.     /*fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);*/
  552.       }
  553.     }
  554. #ifdef MULTIPROC
  555.     }
  556. #endif
  557.     if(!(onRow % 16)) {
  558.       if(outtext)
  559.     fputc('\n',outfile);
  560.       fprintf(stderr,"Done row %d/%d \n",onRow,rows);
  561.     }
  562.   }
  563.  
  564.   if(viewprogram)
  565.     pclose(outfile);
  566.   else if(outname)
  567.     fclose(outfile);
  568.  
  569.   return 0;
  570. }
  571.